Density-Matrix-Power Functional: Performance for Finite Systems and the 

Homogeneous Electron Gas 
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An exchange correlation energy lunctional involving fractional power of the one-body reduced 
density matrix [Phys. Rev. B 78, 201103 (2008)] is applied to finite systems and to the homogeneous 
electron gas. The performance of the functional is assessed for the correlation and atomization 
energies of the molecules contained in the G2 set and for the correlation energy of the homogeneous 
electron gas. High accuracy is found for these two very different types of systems. 



For the past 40 years, density functional theory (DFT) 
developed into one of the most successful theories in the 
study of the electronic structure of atoms, molecules, and 
periodic solids. At its heart lies the exchange-correlation 
(xc) functional, for which many approximations have 
been proposed. The simplest functionals, that depend 
only on the density (the local density approximation - 
LDA), or on the density and its gradients (the generalized 
gradient approximation - GGA), give a very satisfactory 
description of many ground-state properties. However, 
they still fail to reach chemical accuracy for some impor- 
tant quantities like reaction or atomization energies. To 
remedy this situation hybrid functionals were introduced, 
the first and most widely used example being the B3LYP 
functional fiiS This functional is able to reproduce experi- 
mental atomization energies within about 10% error. Al- 
though the atomization energies obtained using B3LYP 
are in good agreement with experiments, the absolute 
correlation energies, an accurate description of which can 
be thought of as a test for the quality of any approximate 
functional, exhibit a sizeable error (up to 400%) 4 This 
is not a surprise since experimentally one normally mea- 
sures energy differences, and it is these quantities that 
functionals like B3LYP are designed to reproduce. Ac- 
curate correlation energies for finite systems can be ob- 
tained by going beyond the DFT framework, for instance 
by using M0ller-Plesset second-order perturbation the- 
ory (MP2) or the coupled cluster method with singles, 
doubles and perturbativc triples [CCSD(T)]. However, 
these methods are computationally too expensive to be 
applied to realistic systems of ever growing complexity: 
bio- molecules, large clusters, nano-devices to name but a 
few examples. 

Recently, reduced density matrix functional theory 
(RDMFT) has appeared as an alternative approach to 
handle complex systems. It has shown great potential for 
improving upon DFT results for finite systems. RDMFT 
uses the one-body reduced density matrix (1-RDM), 7, as 
the basic variable4i^ This quantity, for the ground state. 



is determined through the minimization of the total en- 
ergy functional, under the constraint that 7 is ensemble 
A^-rcprcsentable. The total energy as a functional of 7 
can be expressed as 
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where, p{y) — the electron density — is the diagonal of 
the 1-RDM and v{t) is the external potential. The first 
two terms in Eq. ([1]) are the kinetic and external poten- 
tial energies. The electron-electron interaction can be 
cast in the last two terms, the first being the Coulomb 
repulsion energy and E^c the exchange-correlation (xc) 
energy functional. In principle, Gilbert's^ generalization 
of the Hohcnberg-Kohn theorem to the 1-RDM guaran- 
tees the existence of a functional i?u[7] whose minimum 
yields the exact 7 and the exact ground-state energy 
of the systems characterized by the external potential 
v{y). In practice, however, the xc energy is an unknown 
functional of the 1-RDM and needs to be approximated. 
In the last years, a plethora of approximate function- 
als have been introduced An assess- 
ment of the relative performance of these functionals for 
a large set of atoms and molecules reveals that the so 
called BBC3H and PNOFOM, functionals yield results for 
molecular systems, with an accuracy comparable to the 
MP2 method^'ii^'ii^'ia 

The situation for extended systems is somewhat more 
complicated. When applied to the simplest system, 
the homogeneous electron gas, these functionals lead to 
rather inaccurate correlation energiesi^i Moreover, for 
periodic solids, the aforementioned functionals fail to 
reproduce the fundamental gaps for insulators^ (band 
as well as Mott). A new functional was introduced by 
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Sharma et al. to solve this problcmi^S It reads 

,3 , l7"(r,r')p 



(2) 



where, 7" mdicates the power used in the operator sense 
i.e., 



'(r,r') = ^K-r0*(r') 



(3) 



Here (r) denote the natural orbitals and nj their occu- 
pation numbers. The functional in Eq. ^ was shown to 
perform very well for solidsi^ As it involves the power 
of the density matrix we refer to it as "power func- 
tional" in the following. The power a lies in the interval 
1/2 < a < 1. In the limit a — 1, minimization of the to- 
tal energy functional H]) yields the Hartree-Fock energy, 
while the case a = 1/2 corresponds to the Miiller func- 
tional which tends to over correlate i^s^i Hence the power 
a plays the similar role as the mixing parameter in the 
hybrid functionals of the DFT. 

The situation as it stands is that most of the existing 
xc functionals of RDMFT are designed for finite systems 
and perform quite poorly for solids and the HEG. The 
power-functional, on the other hand, is designed for the 
case of solids, but has not yet been applied to the HEG 
or finite systems. The objective of the present work is 
to fill this gap, in the pursuit of a functional form that 
works equally well for both finite and extended systems. 
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FIG. 1: (Color online) Correlation energy as a function of 
the Wigner-Seitz radius for the homogeneous electron gas. 
RDMFT results are obtained using various approximations 
to the xc functional. Monte Carlo results are taken from 
Ref. m. 

It is difficult to overstate the importance of the HEG 
in the development of many-body theories. It is not only 
used as a benchmark, but also acts as a reference system 
for the design of xc functionals. Within DFT, the LDA 
is perhaps the most successful example of this. Further- 
more, results for the HEG give us valuable indications 
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FIG. 2: (Color online) Percentage deviation of the correlation 
energy, obtained using various 1-RDM functionals, from the 
exact CCSD(T) results. 



on how the theory will perform especially for metallic 
systems. 

With this in mind, we first compare the relative per- 
formance of various RDMFT functionals in reproducing 
the correlation energy of the HEG of various densities. 
Fig. [1] is a plot of the correlation energy versus the den- 
sity parameter rs for a variety of approximate function- 
als compared to exact Monte-Carlo valuesj^ Clearly, the 
Miiller^ CHF— and CGAi2i functionals perform poorly 
over the whole range of r^. The BBGi^ and PNOFii 
functionals are more reasonable but still far from the 
true result. Encouragingly, we find that for a between 
0.55 and 0.58, the power-functional lies very close to the 
Monte Carlo results and possesses a good low density 
limit, making it one of the best 1-RDM functionals for 
the HEG. 

Since the power-functional performs well for the HEG 
at small values of Ts and for periodic solidsj^ it is worth- 
while to investigate its behavior for finite systems. First 
we performed a free optimization of the parameter a us- 
ing a set of 54 molecules and radicals. These molecules 
form a subset of the G2 test set of moleculesi^i^ For 
this optimization two different basis sets are employed, 
namely the 6-31G* and the cc-pVDZ. The optimal value 
of a that minimizes the error in the correlation energy 
for this set of molecules is 0.578. This value is essen- 
tially the same for both kinds of basis set used in the 
present work. It is interesting to note from Fig. [1] that 
the value of a which best reproduces the Monte-Carlo 
data for of interest for most atoms and molecules is 
about 0.55. The optimal value obtained for the set of 
molecules, a = 0.578, is quite close to this. 

Having determined the optimal value for the parame- 
ter a, we compare this approach with different RDMFT 
functionals. Fig. [5] is the plot of relative error in the 
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FIG. 3: (Color online) Atomization energies for the G2 set 
of molecules calculated using the BBC3, PNOFO and power- 
functional vs the CCSD(T) reference atomization energies. 



correlation energy (S) given by 
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where Ec is the correlation energy obtained with the 
method under consideration, Ec°^^ the reference correla- 
tion energy which is obtained with the CCSD(T) method, 
and iVmoi the number of systems included in the calcula- 
tion. It is apparent from Fig. [2]that the power-functional 
is very good in determining the correlation energy of fi- 
nite systems. In fact, we find that this very simple func- 
tional with q;=0.578 (same value of a was used for the full 
G2 set) is of similar quality as the BBC3 or the PNOFO 
which arc to-date the best RDMFT functionals as far as 
finite systems are concernedi^ i^^'^^i^^'^^i^^ 

The accurate determination of atomization energies is 
important for the calculation of the energetics of any 
chemical reaction. Hence this important quantity also 
acts as a test for the quality of an approximate functional. 
In Fig. [3] we plot the atomization energies of molecules 
of the entire G2 set obtained with various approximate 
functionals within RDMFT versus the reference value de- 
termined using the CCSD(T) method. 

It is clear from Refi and Fig. [2] that BBC3, PNOFO 
and the power-functional are the the most accurate 
among the xc functionals we considered, hence in Fig. 
121 we only compare these three functionals. It is clearly 
visible that the power-functional is comparable in accu- 
racy to the BBC3 functional, while PNOFO is slightly 
better than the two. 

The successful prediction of properties of molecules at 
equilibrium does not necessarily imply a good perfor- 
mance for strongly distorted molecules, the dissociation 
limit being one such example. For a stretched molecule, 
not only the total energy has to be equal to the sum of 
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FIG. 4: (Color online) Total energy of the H2 molecule (in 
a.u.) vs the distance (in A) between the two hydrogen atoms. 
RDMFT results are obtained using various approximations to 
the xc functional. For reference the configuration interaction 
results are also presented (black line). 



the energies of the fragments but also the occupations 
of the natural orbitals have to be correct. The simplest 
example is perhaps the II2 molecule. If we take the hy- 
drogen atoms far apart, the total energy should go to -1 
a.u. and the four occupied spin-orbitals made from the 
Hydrogen Is states, have to be occupied by half an elec- 
tron each. Many DFT and RDMFT functionals fail to 
reproduce the correct dissociation of this systemj ^^i^^ 

In Fig.m we plot the H2 dissociation curve obtained us- 
ing various 1-RDM functionals, together with the exact 
curve obtained through a full configuration-interaction 
calculation. The BBC3 functional is designed with the 
dissociation limit in mind and it yields the most accurate 
results for the H2 dissociation. The failure of PNOFO 
and GU functionals is two fold; as the distance between 
the two H atoms increases, the energy deviates strongly 
from the exact value. Second, at a sufficiently large dis- 
tance between the two H atom (6 A) two of the bond- 
ing orbitals still have occupancy 0.86 and the other two 
0.14. Both the Miiller and the power- functional with 
a — 0.578 give the correct occupancy of ~0.5 for all four 
bonding spin-orbitals, but the dissociation energy is un- 
derestimated by the former of these functionals and over- 
estimated by the later. If, on the other hand, the value 
of a is changed to 0.525 the power-functional describes 
the H2 dissociation curve accurately, with an accuracy 
comparable to the BBC3 functional. 

In all the examples studied, it is clear that the sim- 
ple form of the power-functional suffices to obtain very 
good electronic properties. However, we are faced with a 
problem: the optimal value for a varies from one kind of 
system to another. In fact, we find a = 0.65 for solids. 
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a = 0.55 for the HEG, a = 0.578 for molecules at equilib- 
rium, and a ~ 0.525 for stretched H2. Although one can 
use different values of a for different materials, it would 
be desirable to have a unique method to determine the 
system dependent value of a. In this regard, one could 
make a itself a functional of the 1-RDM and optimize it 
in as ab-initio manner for each system. Many other im- 
provements of the power-functional are also conceivable; 
for example it could be a basis for sophisticated correc- 
tions like those of Gristenko et. alM- and/or removal of 
self-interaction terms. How these modifications affect the 
quality of the power-functional will be the subject of fu- 
ture studies. 

In summary, we have used the recently proposed 
power-functional to calculate the correlation energy of 
the HEG, G2 test set of molecules, and stretched H2. For 



the case of molecules, we also determined atomization en- 
ergies. Our results show that the power-functional, orig- 
inally proposed for solids, also performs very well for the 
HEG and finite systems. However, the optimal value of 
a for all three cases is different. The encouraging results 
of the present work point to the fact that this family of 
approximations is an interesting path for the future de- 
velopment of approaches, within RDMFT, to accurately 
describe electronic correlations. 
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gemeinshaft (SPP 1145), and by the EC Network of 
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Coimbra. 
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